%parameter
%clc;clear;
cd /Volumes/Samsung_T5/results_hour_mean/data
r1=[130:300];r2=[305:360];
num1=2;num2=1;num3=744;%time
num4=300-130+1;num5=360-305+1;num6=84;
depth=[1,3,7,12,...
       18,24,34,44,54,64,74,84,94,104,...
       114,124,134,144,154,164,174,184,194,204,...
       214,224,234,244,254,264,274,284,294,304,...
       314,324,334,344,354,364,374,384,394,406,...
       418,432,446,462,478,496,514,534,554,574,...
       594,614,639,664,689,714,744,774,814,854,...
       904,954,1004,1054,1114,1174,1244,1334,1444,1584,...
       1754,1954,2254,2654,3054,3554,4054,4654,5254,5854]; 
%load the corresponding lat/lon; lat_int lon_int
load('lat_int');load('lon_int');
lat_int=lat_int';lon_int=lon_int';
n_x=num4;n_y=num5;
for i=1:n_y
    for j=1:n_x
        if lon_int(i,j)<0
            lon_int(i,j)=lon_int(i,j)+360;
        end
    end
end

%grid distance
count=1;
for i=2:n_x
    x_dis1(:,count)=(lon_int(:,i)-lon_int(:,i-1)).*100000;
    x_dis2(:,count)=(lat_int(:,i)-lat_int(:,i-1)).*100000;
    x_dis(:,count)=sqrt(x_dis1(:,count).*x_dis1(:,count)+x_dis2(:,count).*x_dis2(:,count)); %transfer to m
    count=count+1;
end
x_dis(:,num4)=x_dis(:,num4-1);

x1_dis(1:num5,1)=1;
for i=2:n_x
    x1_dis(:,i)=x_dis(:,i-1)+x1_dis(:,i-1);
end

count=1;
for i=2:n_y
    y_dis1(count,:)=(lon_int(i,:)-lon_int(i-1,:)).*100000;
    y_dis2(count,:)=(lat_int(i,:)-lat_int(i-1,:)).*100000;
    y_dis(count,:)=sqrt(y_dis1(count,:).*y_dis1(count,:)+y_dis2(count,:).*y_dis2(count,:)); %transfer to m
    count=count+1;
end
y_dis(num5,:)=y_dis(num5-1,:);

y1_dis(1,1:num4)=1;
for i=2:n_y
    y1_dis(i,:)=y_dis(i-1,:)+y1_dis(i-1,:);
end

for i=1:83
    zz(i)=depth(1,i+1)-depth(1,i);
end

for i=1:84
    XXX(:,:,i)=x1_dis(:,:);
    YYY(:,:,i)=y1_dis(:,:);
    ZZZ_new(1:n_y,1:n_x,85-i)=depth(1,i);
end

%choose the section along 53N
fid7=fopen('nakano_topo_rnewb2.grd','r');
topo=fread(fid7,800*800,'real*4');
topo_1=reshape(topo,800,800);
topo_1(topo_1<-10000)=nan;
topo_1=topo_1*10;
topo2(:,:)=topo_1(r1,r2)';
fclose(fid7);

%%cross- and along- unit vector
[gx,gy]=gradient(topo2);
gx=gx./x_dis;gy=gy./y_dis;
length_n=sqrt(gx.*gx+gy.*gy);
gx=gx./length_n;gy=gy./length_n; %gradient direction small--->large; right%downward:positive gradient
ay=sqrt(1./((gy.*gy)./(gx.*gx)+1));%tangent direction
ax=(-1.*ay.*gy)./gx;
[fig_xx,fig_yy]=meshgrid(1:171,1:56);

clear xxx_p yyy_p ZZZ zzz_p
xxx_p=fig_xx(:);yyy_p=fig_yy(:);
ZZZ=gx(:,:);
zzz_p=ZZZ(:);
xxx_p(isnan(ZZZ))=[];
yyy_p(isnan(ZZZ))=[];
zzz_p(isnan(ZZZ))=[];
gx_new(:,:)=griddata(xxx_p,yyy_p,zzz_p,fig_xx,fig_yy);

clear xxx_p yyy_p ZZZ zzz_p
xxx_p=fig_xx(:);yyy_p=fig_yy(:);
ZZZ=gy(:,:);
zzz_p=ZZZ(:);
xxx_p(isnan(ZZZ))=[];
yyy_p(isnan(ZZZ))=[];
zzz_p(isnan(ZZZ))=[];
gy_new(:,:)=griddata(xxx_p,yyy_p,zzz_p,fig_xx,fig_yy);

clear xxx_p yyy_p ZZZ zzz_p
xxx_p=fig_xx(:);yyy_p=fig_yy(:);
ZZZ=ax(:,:);
zzz_p=ZZZ(:);
xxx_p(isnan(ZZZ))=[];
yyy_p(isnan(ZZZ))=[];
zzz_p(isnan(ZZZ))=[];
ax_new(:,:)=griddata(xxx_p,yyy_p,zzz_p,fig_xx,fig_yy);

clear xxx_p yyy_p ZZZ zzz_p
xxx_p=fig_xx(:);yyy_p=fig_yy(:);
ZZZ=ay(:,:);
zzz_p=ZZZ(:);
xxx_p(isnan(ZZZ))=[];
yyy_p(isnan(ZZZ))=[];
zzz_p(isnan(ZZZ))=[];
ay_new(:,:)=griddata(xxx_p,yyy_p,zzz_p,fig_xx,fig_yy);
%%
%figure(1)
%contour(LON,LAT,c2,'linewidth',0.3,'color','k');hold on;quiver(LON,LAT,gx,gy);%verify the normal direction
clear test1 test2 test3 test4
pick_lat=53.5;
figure(2);pcolor(topo2);shading interp;
colormap(flipud(othercolor('Accent5')));
hold on;[cc1,cch1]=contour(lat_int,[48 49 50 51 52 53 54 55],'showtext','on','linewidth',3,'color','r','linestyle',':');
clabel(cc1,cch1,'fontsize',10,'labelspacing',500,'FontName','Arial Black','color','k');
hold on;[cc2,cch2]=contour(lon_int,[141 142 143 144 145 146 147 148],'showtext','on','linewidth',3,'color','r','linestyle',':');
clabel(cc2,cch2,'fontsize',10,'labelspacing',500,'FontName','Arial Black','color','k');
hold on;[cc3,cch3]=contour(topo2,[100 150 200 250],'showtext','on','linewidth',6,'color','r');
clabel(cc3,cch3,'fontsize',30,'labelspacing',500,'FontName','Arial Black','color','k');
[test1,test2]=find(abs(lat_int-pick_lat)<0.01);
num_size=size(test1,1);
% figure(2);test3(:,:)=u_rho_jan(test1,test2,:);
for j=1:num_size-1
    if test1(j,1)==test1(j+1,1)
        test1(j,1)=nan;
        test2(j,1)=nan;
    end
end
count=1;
for j=1:num_size
    if ~isnan(test1(j,1))
        test3(count,1)=test1(j,1);
        test4(count,1)=test2(j,1);
        count=count+1;
    end
end
num_size1=size(test4,1);
figure(2);hold on;scatter(test4(:,1),test3(:,1),200,'filled','r');
set(gca,'xticklabel',[],'yticklabel',[],'linewidth',8);
set(gcf,'unit','centimeters','position',[0,0,40,40]);
%print('location2','-dtiff','-r600');
close(figure(2));

%%
% %offshore distance
for i=1:num_size1
    lat_lat(i,1)=lat_int(test3(i,1),test4(i,1),1);
    lon_lat(i,1)=lon_int(test3(i,1),test4(i,1),1);
end

off_dis(1,1)=0;
count=2;
%for j=35:num_size1 %53N

for j=39:num_size1
    delta_lon=lon_lat(j,1)-lon_lat(j-1,1);
    dis=cos(lat_lat(j,1)/180*pi)*cos(lat_lat(j-1,1)/180*pi)*(1-cos(delta_lon/180*pi))/2;
    off_dis(count,1)=6371*2*asin(sqrt(dis));
    count=count+1;
end
% for j=39:num_size1 %53.5N
%     off_dis1(count,1)=(lon_lat(j,1)-lon_lat(j-1,1)).*100000;
%     off_dis2(count,1)=(lat_lat(j,1)-lat_lat(j-1,1)).*100000;
%     off_dis(count,1)=sqrt(off_dis1(count,1).*off_dis1(count,1)+off_dis2(count,1).*off_dis2(count,1)); %transfer to m
%     count=count+1;
% end

size1=size(off_dis);
off1_dis(1,1)=0;
for i=2:size1+1
    off1_dis(i,1)=off_dis(i-1,1)+off1_dis(i-1,1);
end
%off1_dis=off1_dis./1000;
%off1_dis=roundn(off1_dis,-1);
cd /Volumes/Samsung_T5

